Digitalisierung der Trainingsdaten

#---- Digitalisierung der Trainingsdaten ----
# Muß natürlich nur einmal gemacht werden
if (digitize) {
  source(file.path(root_folder, "src/functions/digitize.R"))
} else {
  train_areas_2019 = readRDS(paste0(envrmt$path_data,"train_areas_2019.rds"))
  train_areas_2020 = readRDS(paste0(envrmt$path_data,"train_areas_2020.rds"))
}

Überwachte mit Klassifikation Random Forest

## ---- Überwachte  mit Klassifikation Random Forest ----
## Hier wird der Random Forest über das Utility Paket caret aufgerufen

# Setzen eines "seed" ermöglicht reproduzierbaren Zufall
set.seed(123)

# Zufälliges Ziehen von 25% der Daten
trainDat_2019 =  train_areas_2019[createDataPartition(train_areas_2019$class,list = FALSE,p = 0.25),]
trainDat_2020 =  train_areas_2020[createDataPartition(train_areas_2020$class,list = FALSE,p = 0.25),]
#  Response-Variable (=Spalte "class") muss den Datentyp "factor" haben
trainDat_2019$class <- as.factor(trainDat_2019$class)
trainDat_2020$class <- as.factor(trainDat_2020$class)

# Einstellungen Modelltraining: cross-validation, 10 Wiederholungen
ctrlh = trainControl(method = "cv",
                     number = 10,
                     savePredictions = TRUE)

#--- random forest model training
cv_model_2019 = train(trainDat_2019[,2:11], # in den Spalten 2 bis 20 stehen die Trainingsdaten (Prediktoren genannt)
                      trainDat_2019[,1],         # in der Spalte 1 steht die zu klassizierende Variable (Response genannt)
                      method = "rf",             # Methode hier rf für random forest
                      metric = "Kappa",          # Qualitäts/Performanzmaß KAppa
                      trControl = ctrlh,         # obig erzeugte Trainingssteuerung soll eingelsen werden
                      importance = TRUE)         # Die Bedeung der Variablen wird mit abgespeichert
# Klassifikation (auch Vorhersage genannt)
prediction_rf_2019  = raster::predict(pred_stack_2019 ,cv_model_2019, progress = "text")

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================================| 100%
## 

#--- 2020
# Einstellungen wie 2019
cv_model_2020 = train(trainDat_2020[,2:11],
                      trainDat_2020[,1],
                      method = "rf",
                      metric = "Kappa",
                      trControl = ctrlh,
                      importance = TRUE)
prediction_rf_2020  = raster::predict(pred_stack_2020 ,cv_model_2020, progress = "text")

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================================| 100%
## 

saveRDS(prediction_rf_2019,paste0(envrmt$path_data,"pred_rf_2019.rds"))
saveRDS(prediction_rf_2020,paste0(envrmt$path_data,"pred_rf_2020.rds"))
prediction_rf_2019 = readRDS(paste0(envrmt$path_data,"pred_rf_2019.rds"))
prediction_rf_2020 = readRDS(paste0(envrmt$path_data,"pred_rf_2020.rds"))

#---- Visualisierung mit tmap qtm() ----

# (interaktiv = "view", statisch = "plot")
tmap_mode("plot")

## tmap mode set to plotting

w1 = qtm(prediction_rf_2019,projection = 32632,inner.margins=0.01) +
  tm_legend(scale= 0.5,
            legend.outside=T,
            title = "RF 2019",
            title.size = 1.0) +
  tm_grid()

w2 = qtm(prediction_kmeans_2019$map,projection = 32632,inner.margins=0.01) +
  tm_legend(scale= 0.5,
            legend.outside=T,
            title = "kmeans 2019",
            title.size = 1.0) +
  tm_grid()

w3 = qtm(prediction_rf_2020,projection = 32632,inner.margins=0.01) +
  tm_legend(scale= 0.5,
            legend.outside=T,
            title = "RF 2020",
            title.size = 1.0) +
  tm_grid()

w4 = qtm(prediction_kmeans_2020$map,projection = 32632,inner.margins=0.01) +
  tm_legend(scale= 0.5,
            legend.outside=T,
            title = "kmeans 2020",
            title.size = 1.0) +
  tm_grid()
# ohne jede Anordnung
w2;w4;w1;w3

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

# mit der bordeigenen tmap fnktion tmap_arrrange()
tmap_arrange(w1, w3, w2, w4,asp = NA)

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

# mit dem paket gid (universell nutzbar um graphische ausgabe zu gestalten)
# hier die aufteilung wi tmap_arrange()
library(grid)

## 
## Attaching package: 'grid'

## The following objects are masked from 'package:terra':
## 
##     arrow, depth

grid.newpage()
print(w2, vp=viewport(0.25, 0.25, width = 0.5, height = 0.5))

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

print(w1, vp=viewport(0.25, 0.75, width = 0.5, height = 0.5))

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

print(w4, vp=viewport(0.75, 0.25, width = 0.5, height = 0.5))

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

print(w3, vp=viewport(0.75, 0.75, width = 0.5, height = 0.5))

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

# mit dem paket cowplot
library(cowplot)
plot_grid(tmap_grob(w1),tmap_grob(w3),tmap_grob(w2),tmap_grob(w4),
          labels = c('FIG A', 'FIG B', 'FIG C', 'FIG D'),
          label_size= 10, vjust = 5.)

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

## Note that the aspect ratio for which the grob has been generated is 1.4

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

## Note that the aspect ratio for which the grob has been generated is 1.4
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

## Note that the aspect ratio for which the grob has been generated is 1.4
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

## Note that the aspect ratio for which the grob has been generated is 1.4

### Maximum Likelihood Classification

# ---- Maximum Likelihood Classification ----
# superClass() Funktion aus dem Paket RSToolbox umwandeln der Tabelle in das
# geforderte SpatialdataPoint Objekt Identifikation der Koordinaten
sp::coordinates(trainDat_2019) = ~X+Y
sp::coordinates(trainDat_2020) = ~X+Y
# Zuweisen der korrekten Projektion (hier aus dem zugehörigen Raster-Stack)
crs(trainDat_2019) = crs(pred_stack_2019)
crs(trainDat_2020) = crs(pred_stack_2020)

# superClass method "mlc" berechnet die Statisik und klassifiziert in einem aufruf
prediction_mlc_2019       <- superClass(pred_stack_2019, trainData = trainDat_2019[,1:11], responseCol = "class",
                                        model = "mlc", tuneLength = 1, trainPartition = 0.5)

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

prediction_mlc_2020       <- superClass(pred_stack_2020, trainData = trainDat_2020[,1:11], responseCol = "class",
                                        model = "mlc", tuneLength = 1, trainPartition = 0.5)

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.

saveRDS(prediction_mlc_2019,paste0(envrmt$path_data,"pred_mlc_2019.rds"))
saveRDS(prediction_mlc_2020,paste0(envrmt$path_data,"pred_mlc_2020.rds"))
prediction_mcl_2019 = readRDS(paste0(envrmt$path_data,"pred_mlc_2019.rds"))
prediction_mlc_2020 = readRDS(paste0(envrmt$path_data,"pred_mlc_2020.rds"))

# Ergebnisse der rf Klassifikation
prediction_mlc_2019

## superClass results
## ************ Validation **************
## $validation
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction clearcut other
##   clearcut       37    11
##   other           5   670
##                                           
##                Accuracy : 0.9779          
##                  95% CI : (0.9643, 0.9873)
##     No Information Rate : 0.9419          
##     P-Value [Acc > NIR] : 2.471e-06       
##                                           
##                   Kappa : 0.8105          
##                                           
##  Mcnemar's Test P-Value : 0.2113          
##                                           
##             Sensitivity : 0.88095         
##             Specificity : 0.98385         
##          Pos Pred Value : 0.77083         
##          Neg Pred Value : 0.99259         
##              Prevalence : 0.05809         
##          Detection Rate : 0.05118         
##    Detection Prevalence : 0.06639         
##       Balanced Accuracy : 0.93240         
##                                           
##        'Positive' Class : clearcut        
##                                           
## 
## *************** Map ******************
## $map
## class      : RasterLayer 
## dimensions : 1130, 1499, 1693870  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 582400, 597390, 5744870, 5756170  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : class 
## values     : 1, 2  (min, max)
## attributes :
##  ID    value
##   1 clearcut
##   2    other

prediction_mlc_2020

## superClass results
## ************ Validation **************
## $validation
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction clearcut other
##   clearcut      114    28
##   other           7   766
##                                           
##                Accuracy : 0.9617          
##                  95% CI : (0.9472, 0.9732)
##     No Information Rate : 0.8678          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.8448          
##                                           
##  Mcnemar's Test P-Value : 0.0007232       
##                                           
##             Sensitivity : 0.9421          
##             Specificity : 0.9647          
##          Pos Pred Value : 0.8028          
##          Neg Pred Value : 0.9909          
##              Prevalence : 0.1322          
##          Detection Rate : 0.1246          
##    Detection Prevalence : 0.1552          
##       Balanced Accuracy : 0.9534          
##                                           
##        'Positive' Class : clearcut        
##                                           
## 
## *************** Map ******************
## $map
## class      : RasterLayer 
## dimensions : 1130, 1499, 1693870  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 582400, 597390, 5744870, 5756170  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : class 
## values     : 1, 2  (min, max)
## attributes :
##  ID    value
##   1 clearcut
##   2    other

## ---- Visualisierung mit mapview ----
mapview::viewRGB(mask*pred_stack_2020, r = 4, g =5, b = 6,maxpixels =  1693870)+
  mapview(mask*prediction_rf_2019 , alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = TRUE) +
  mapview(mask*prediction_rf_2020, alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE) +
  mapview(mask*prediction_mlc_2019$map,alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE) +
  mapview(mask*prediction_mlc_2020$map,alpha.regions = 0.5, maxpixels =  1693870,
          col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE)

## Warning in raster::projectRaster(x, raster::projectExtent(x, crs =
## sp::CRS(epsg3857)), : input and ouput crs are the same

Change Detection Auswertung

## ---- Change Detection Auswertung ----

# Differenzbild random forest
tmap::qtm(prediction_rf_2020 - prediction_rf_2019,title = "kmeans 2019")

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

# ---- CDA Extraktion der Klassennamen ----
categories = levels(prediction_rf_2019)[[1]]$value
categories

## [1] "clearcut" "other"

# ---- Berechnung der Vierfeld Tabelle mit raster Basisfunktion ----
ct = raster::crosstab(prediction_rf_2019,prediction_rf_2020)
rownames(ct) = categories
colnames(ct) = categories
ct %>%
  kbl(caption = "Crosstab 2019 vs 2020",) %>%
  kable_classic(full_width = F)

# add_header_above(c("layer 1 " = 1, "layer 2" = 2))

# ---- kappa ----
# Vergleich der Übereinstimmung unterschiedlicher Klassifikationen (hier MaxLike
# und RF) mit Hilfe diverser Kappa Werte
# https://giswerk.org/doku.php?id=r:r-tutorials:calculatekappa
# 2019
kappa_2019 = kstat(prediction_mlc_2019$map,prediction_rf_2019, perCategory = FALSE)
kappa_2019

##                 K      Kloc    Khisto CA       QA       AA        AD       QD
## overall 0.5285764 0.9298423 0.5684582 50 46.68416 1.752676 0.1322416 1.430924

kappa_2019  %>%
  kbl(caption = "Kappa 2019 MLC vs RF",) %>%
  kable_classic(full_width = F)

# 2020
kappa_2020 = kstat(prediction_mlc_2020$map,prediction_rf_2020,perCategory = FALSE)
kappa_2020

##                 K      Kloc    Khisto CA      QA       AA       AD       QD
## overall 0.6513909 0.8387926 0.7765815 50 43.2518 4.395713 0.844811 1.507672

kappa_2020  %>%
  kbl(caption = "Kappa 2019 MLC vs RF",) %>%
  kable_classic(full_width = F) %>%
  column_spec(1,color = "red", link = "https://giswerk.org/doku.php?id=r:r-tutorials:calculatekappa")

# Erzeugen aller Vergleichsraster der Kontingenztabelle
# Die lapply funktionen sind integrierte FOR Schleifen die über die Liste der
# Kategorien die Funktion changefrom() für die Kreuztabelle anwenden
r = lapply(1:length(categories),
           function(i){lapply(1:length(categories),
                              function(j){changefrom(prediction_rf_2019, prediction_rf_2020, i,j)})})
r

## [[1]]
## [[1]][[1]]
## class      : RasterLayer 
## dimensions : 1130, 1499, 1693870  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 582400, 597390, 5744870, 5756170  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)
## 
## 
## [[1]][[2]]
## class      : RasterLayer 
## dimensions : 1130, 1499, 1693870  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 582400, 597390, 5744870, 5756170  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)
## 
## 
## 
## [[2]]
## [[2]][[1]]
## class      : RasterLayer 
## dimensions : 1130, 1499, 1693870  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 582400, 597390, 5744870, 5756170  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)
## 
## 
## [[2]][[2]]
## class      : RasterLayer 
## dimensions : 1130, 1499, 1693870  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 582400, 597390, 5744870, 5756170  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)

# ---- Visualsierung der Kreuztabellierten Von-Zu-Raster ----
# Plotten der Raster hierzu werden erneut alle Kategorien einzeln geplottet i
# und j sind hilfsvariablen um die korrekten Raster Layer ansprechen zu können.
# t ist eine Hilfvariable um eine Liste für die Ergebnisbilder hochzählen zu
# können
tmap_mode("view")

## tmap mode set to interactive viewing

t=i=j=1 # setze zählvariablen auf 1
m=list() # erzeuge leere Liste für die Karten
for(cat1 in categories)  { # für jede Kategorie
  for(cat2 in categories)  { # mit jeder Kategorie
    # plotte die Karte
    m[[t]]  = tm_shape(st_as_stars(r[[i]][[j]])) +
      tm_raster(col = "layer",palette = "viridis",style = "cat",
                title = if(cat1==cat2) {
                  paste("no changes " , unique(cat1,cat2))
                }
                else if (cat1!=cat2) {
                  paste(cat1," -> ",cat2)
                })
    # zähle die innere schleife hoch
    j = j + 1
    # zähle die ergebnisliste hoch
    t = t + 1
  }
  # innere schleife abgearbeitet setze sie auf  den Anfang zurück
  j = 1
  # zähle die äußere Schleife hoch
  i = i + 1
}
# Interaktive und synchronisierte Karten
tmap::tmap_arrange(m,sync = TRUE)

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)

# ---- Analyse von Veränderungen mit dem paket motif ----
# lsp_compare vergleicht zwei (oder mehr) kategoriale Karten (change detection
# klassifikationen etc.) miteinander und nutzt für die Ausgabe der
# Wahrscheinlichkeiten verschiedener räumlicher Aggregierungsstufen und
# Merkmalsraum-Distanzen um Veränderungswahrscheinlichkeiten zu ermitteln
mrf_compare_2020_2019 = lsp_compare(st_as_stars(mask*prediction_rf_2019),
                                    st_as_stars(mask*prediction_rf_2020),
                                    type = "cove",
                                    dist_fun = "jensen-shannon",
                                    window = win_size,
                                    threshold = 0.9)

## Metric: 'jensen-shannon' using unit: 'log2'.

# Visualisierung der Gesamtwahrscheinlichkeiten Achtung logarithmische Skala

tmap_mode("plot")

## tmap mode set to plotting

tm_compare2 = tm_shape(mrf_compare_2020_2019) +
  tm_raster("dist",
            breaks = c(0, 0.001, 0.01, 0.1, 1.01),
            legend.is.portrait = F,
            palette = "viridis",
            title =  "Distance (JSD)") +
  tm_layout(legend.show = TRUE,
            legend.text.size = 0.3,
            legend.outside = TRUE) +
  tm_legend(scale= 0.5,
            legend.outside=T,
            title = paste("RF 2019 vs 2020 ",win_size*10,"m**2" ),
            title.size = 1.0) +
  tm_grid()

tm_compare2

# ---- Detail-Analyse Teil 1 ----
# Identifikation der Gebiete (Referenz ist win_size) die maximale Veränderungen
# aufweisen im Beispiel soll  "dist" soll größer 0.001 sein (logarithmische
# Skalierung!)
lc_am_compare_sel = st_as_sf(mrf_compare_2020_2019) %>%
  subset(dist > 0.01)

# Sortierung nach Größe
lc_am_compare_sel = lc_am_compare_sel[order(lc_am_compare_sel$dist,decreasing = TRUE), ]
lc_am_compare_sel

## Simple feature collection with 552 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 582400 ymin: 5744920 xmax: 597400 ymax: 5755670
## Projected CRS: WGS 84 / UTM zone 32N
## First 10 features:
##        id na_prop_x na_prop_y      dist                       geometry
## 2107 2107         0         0 0.6196571 POLYGON ((583900 5747420, 5...
## 2291 2291         0         0 0.5594878 POLYGON ((584900 5746670, 5...
## 2290 2290         0         0 0.5212548 POLYGON ((584650 5746670, 5...
## 2470 2470         0         0 0.3918801 POLYGON ((584650 5745920, 5...
## 2471 2471         0         0 0.3341960 POLYGON ((584900 5745920, 5...
## 2650 2650         0         0 0.3203009 POLYGON ((584650 5745170, 5...
## 2231 2231         0         0 0.3052788 POLYGON ((584900 5746920, 5...
## 2531 2531         0         0 0.2897703 POLYGON ((584900 5745670, 5...
## 2116 2116         0         0 0.2495521 POLYGON ((586150 5747420, 5...
## 2649 2649         0         0 0.2460312 POLYGON ((584400 5745170, 5...

##- Visualsierung
# Plotten der top ten gebiete wir nutzen die selbst geschriebene Funktion tm_plot()
# siehe src/functions/fun_tmap.R
tm_plot(sel=lc_am_compare_sel[1:nrow(lc_am_compare_sel),],ov = T)

## tmap mode set to plotting

## Legend labels were too wide. Therefore, legend.text.size has been set to 0.4. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

# als interaktive karte
tm_plot(sel=lc_am_compare_sel[1:nrow(lc_am_compare_sel),],vt = "view",ov = T)

## tmap mode set to interactive viewing

# ---- Visualsierng Detail-Analyse ----
# Extraktion nach der sortierten Liste lc_comapare_sel  werden die top 10 change
# detection hotspots Daten extrahiert und in die liste lc geschrieben
lc = lapply(seq(1:10),function(i){
              lsp_extract(c(st_as_stars(prediction_rf_2019),
                            st_as_stars(prediction_rf_2020)),
                          window = win_size,
                          id = lc_am_compare_sel$id[i])
            })

# Erzeugen der top 3 Change Detection Hotspots Karten-Ausschnitte
tm_lc = lapply(seq(1:3),function(i){
               tm_plot(lc[[i]])
             })

## tmap mode set to plotting

## tmap mode set to plotting
## tmap mode set to plotting

# plotten der erzeugten Karten
tmap_arrange(tm_lc)

#####